/****************************************************************************
 * filterspecial.cpp
 * C++ source file generated from polya project.
 *
 * Release Date: 02/20/2013
 * Author: Min Dong
 *
 * Module version: 0.1.022212
 ****************************************************************************/


#include <math.h>
#include <vector>
#include <map>
#include <iostream>
#include <fstream>
#include<sstream>
#include<stdlib.h>
#include <string.h>
using namespace std;

#define POOL_SIZE 1024
struct parameter
{
	string name;
	string value;
};
struct position
{
	int start;
	int end;
};


/* info struct is used to store additional infomation for FASTA and FASTQ format,
 * Example : FASTQ format
 *  @HWI-ST741:203:D0P2VACXX:1:1101:1167:1791 1:N:0:      						  -->  title
 *	NGGGATTTTTTTTTTTTTTTTTTTTTTTTTAAAATACTTTCGCTACAAATAAAATTTTTGGTCGCATAAATAACC
 *	+ 																		      --> plus_head
 *	#1:BDDDDFDDDDIIIIIIID@68?##################################################   --> plus_content
 *	In case FASTA format, plus_head and plus_content are not used.
 */
struct info
{
	char * title;
	char * plus_head;
	char * plus_content;
	info(){
		title = NULL;
		plus_head = NULL;
		plus_content = NULL;
	}
	void init(int max_size){
		title = new char[max_size];
		plus_head = new char[max_size];
		plus_content = new char[max_size];
	}
	void free() {
		if (title) delete title;
		if (plus_head) delete plus_head;
		if (plus_content) delete plus_content;
	}
};
#define TYPE_NONE 0
#define TYPE_0TR_1TL 1   //tr_num==0 && tl_num==1
#define TYPE_1TR_0TL 2   //tl_num==0 && tr_num==1
#define TYPE_1TR_1TL 3   //tl_num==1 && tr_num==1

#define BLOCK 64
__global__ void kernel_findp(char * lines, int lines_count , int line_lenght, int * type,int * lens, char ptag ,int mismatch,
		int scangap ,int tail_length, int head_gap, float * true_identity, position * pos_list, int * pos_count);
__global__ void kernel_finds(char * lines, int lines_count , int line_lenght, int * type,int * lens, char ptag ,int mismatch,
		int scangap ,int tail_length, int head_gap, float * true_identity, position * pos_list, int * pos_count);


int main(int argc, char* argv[])
{
	int default_sign=0;
	int detail=0;
	if(argc<3 || argv[1]=="-h")
	{
		cout<<">***FindTail help information***"<<endl;
		cout<<">for example:"<<endl;
		cout<<">command:findtail --input_file a.fastq --seqlength 500 --endgap 10 --taillength 15 --identity 80 --ptype T --stype A --output_format fasta --output_type n --d > outputfile"<<endl;
		cout<<">***Version information***"<<endl;
		cout<<"FindTail Version: 0.11.02202013"<<endl;
		cout<<"Author: Min Dong"<<endl;
		cout<<"release at 02/20/2013"<<endl;
		cout<<"***********************************help***************************************"<<endl;
		cout<<">Search module version:0.27.0101111"<<endl;
		cout<<">--input_file: input file name, either in FASTA or FASTQ format;"<<endl;
		cout<<">--seqlength: the maximum sequence length in the input file;"<<endl;
		cout<<">--endgap: the minimum distance between the sequence front end and the tail start or between the sequence back end and the tail end;"<<endl;
		cout<<">--taillength: the minimum length of a tail for detection;"<<endl;
		//cout<<">--mismatch :allow mismatch in tail sequence,0 means perfect match;"<<endl;
		cout<<">--identity: the minimum identity percentage [between 50 and 100] required for a valid tail, where 100 means perfect match;"<<endl;
		cout<<">--ptype: the primary tail type (e.g., poly(A)/(T)/(G)/(C)) for your search;"<<endl;
		cout<<">--stype: the secondary tail type (e.g., poly(T)/(A)/(C)/(G)) for your search;"<<endl;
		cout<<">--output_format: the output format, either in FASTA or FASTQ;"<<endl;
		cout<<">--output_type: the output type that you can choose: [n|pl|pr|sl|sr|pd|pd|sd|np|ns|none]"<<endl;
		cout<<"           n:only output statictic result;"<<endl;
		cout<<"           pl: output primary tail in left region (the distance is defined in --endgap);"<<endl;
		cout<<"           pr: output primary tail in right region (the distance is defined in --endgap);"<<endl;
		cout<<"           sl: output secondary tail in left region (the distance is defined in --endgap);"<<endl;
		cout<<"           sr: output secondary tail in right region (the distance is defined in --endgap);"<<endl;
		cout<<"           pd: output primary tail in left region and right region (the distance is defined in --endgap);"<<endl;
		cout<<"           sd: output secondary tail in left region and right region (the distance is defined in --endgap);"<<endl;
		cout<<"           np: output sequence reads without the primary tail detected;"<<endl;
		cout<<"           ns: output sequence reads without the secondary tail detected;"<<endl;
		cout<<"           none: output sequence reads without both primary and secondary tails detected;"<<endl;
		cout<<">--d: the last parameter,if it exist,the output information contains motif length and identity value;"<<endl;
		cout<<"******************************************************************************"<<endl;

		return 0;
	}

	if(argc >1 && argc<19 && argv[1]!="-h")
	{
		default_sign=1;
		cout<<"please input correct parameters."<<endl;
		return 0;
	}
	vector<parameter> p_container;
	string input_filename;
	int dbfsequenceboundary = 500;
	int tail_length = 15;
	int head_gap = 10;
	int mismatch =2;

	string ptag = "T";
	string stag = "A";
	string output_s = "n";
	string output_format = "fasta";
	if(default_sign==0)
	{
		parameter p;
		//for(int i=1;i<19;i=i+2)
		//{
		//    p.name = argv[i];
		//    p.value=argv[i+1];
		//    //cout<<p.name<<" "<<p.value<<endl;
		//    p_container.push_back(p);

		//}
		////set detail value
		//if(argc==20)
		//{
		//    if(string(argv[19])=="--d")
		//    {
		//        detail=1;
		//    }
		//
		//}
		if(argc==20)
		{
			for(int i=1;i<20;i++)
			{
				if(argv[i][0]=='-')
				{
					if(argv[i][2]!='d')
					{
						p.name = argv[i];
						p.value=argv[i+1];
					}
					else
					{
						p.name = argv[i];
						p.value="none";
					}

				}

				//cout<<p.name<<" "<<p.value<<endl;
				p_container.push_back(p);

			}
		}
		else if(argc==19)
		{
			for(int i=1;i<19;i=i+2)
			{
				p.name = argv[i];
				p.value=argv[i+1];
				//cout<<p.name<<" "<<p.value<<endl;
				p_container.push_back(p);

			}
		}


		//set configuration
		for(int j=0;j<p_container.size();j++)
		{
			if(p_container[j].name=="--d")
			{
				detail=1;
			}
			if(p_container[j].name=="--input_file")
			{
				input_filename =p_container[j].value;
			}
			if(p_container[j].name=="--seqlength")
			{
				dbfsequenceboundary =atoi(p_container[j].value.c_str());
			}
			if(p_container[j].name=="--endgap")
			{
				head_gap =atoi(p_container[j].value.c_str());
			}
			if(p_container[j].name=="--taillength")
			{
				tail_length =atoi(p_container[j].value.c_str());
			}
			/*if(p_container[j].name=="--mismatch")
            {
                mismatch =atoi(p_container[j].value.c_str());
            }*/
			if(p_container[j].name=="--identity")
			{

				float temp = atoi(p_container[j].value.c_str());
				mismatch = (100 - temp)/100 * tail_length;

				if(mismatch==0)
				{
					mismatch=1;
				}
				if(mismatch<2)
				{
					mismatch=2;
				}
				if(mismatch>2)
				{
					mismatch=2;
				}
				//cout<<mismatch<<endl;
			}
			/*if(p_container[j].name=="--tailtype")
            {
                tag =p_container[j].value;
            }*/
			if(p_container[j].name=="--ptype")
			{
				ptag =p_container[j].value;
			}
			if(p_container[j].name=="--stype")
			{
				stag =p_container[j].value;
			}
			if(p_container[j].name=="--output_format")
			{
				output_format =p_container[j].value;

			}
			if(p_container[j].name=="--output_type")
			{
				output_s = p_container[j].value;
				if(p_container[j].value == "pl")
				{
					output_s ="tl";
				}
				if(p_container[j].value == "pr")
				{
					output_s ="tr";
				}
				if(p_container[j].value == "sl")
				{
					output_s ="al";
				}
				if(p_container[j].value == "sr")
				{
					output_s ="ar";
				}
				if(p_container[j].value == "pd")
				{
					output_s ="td";
				}
				if(p_container[j].value == "sd")
				{
					output_s ="ad";
				}
				if(p_container[j].value == "np")
				{
					output_s ="nt";
				}
				if(p_container[j].value == "ns")
				{
					output_s ="na";
				}

			}

		}

		int scangap= 2;
		int div = 10;


		ifstream infile(input_filename.c_str());
		if (!infile)
		{
			cout<<"can't find the polya site file\n";
			return -1;
		}
		else
		{
			char *line = new char[dbfsequenceboundary];
			// Determine format type of input data.
			bool isFastQ = false;
			infile.getline(line,dbfsequenceboundary);
			if (line[0] == '@')
			{
				isFastQ = true;
				div = 4;
			}else {
				isFastQ = false;
				div = 2;
				//fasta must output fasta
				output_format = "fasta";
			}


			info * info_list;
			info_list = new info[POOL_SIZE];
			for (int i = 0 ; i < POOL_SIZE ; i++)
			{
				info_list[i].init(dbfsequenceboundary)	;
			}
			char * lines_p  = new char[dbfsequenceboundary * POOL_SIZE];
			char * lines_s  = new char[dbfsequenceboundary * POOL_SIZE];
			char * temp_title = new char [dbfsequenceboundary];
			int * lens_p= new int[POOL_SIZE];
			int * lens_s= new int[POOL_SIZE];
			float * true_identity_p = new float[POOL_SIZE * 2];
			float * true_identity_s = new float[POOL_SIZE * 2];
			int * type_p = new int[POOL_SIZE];
			int * type_s = new int[POOL_SIZE];
			position * pos_list_p = new position[POOL_SIZE * 2];
			position * pos_list_s = new position[POOL_SIZE * 2];
			int * pos_count_p = new int [POOL_SIZE];
			int * pos_count_s = new int [POOL_SIZE];

			//CUDA declaration
			int * d_lens;
			char * d_lines;
			int * d_type;
			float * d_true_identity;
			position * d_pos_list;
			int * d_pos_count;
			cudaMalloc((void **)& d_pos_count, POOL_SIZE * sizeof(int));
			cudaMalloc((void **)& d_type, POOL_SIZE * sizeof(int));
			cudaMalloc((void **)& d_lens, POOL_SIZE * sizeof(int));
			cudaMalloc((void **) & d_lines,dbfsequenceboundary * POOL_SIZE * sizeof(char));

			cudaMalloc((void **) & d_true_identity, 2 * POOL_SIZE * sizeof(float));
			cudaMalloc((void **) & d_pos_list,  2 * POOL_SIZE * sizeof(position));
			long count =0;

			double t_left =0;
			double t_right=0;
			double a_left=0;
			double a_right=0;
			double t_double=0;
			double a_double=0;
			double tl_ar=0;
			double al_tr=0;
			//infomation for sequence
			//       string title_temp;
			//       string quality_temp;
			//       //quality sign
			//       int quality_sign=0;

			while(!infile.eof())
			{
				int pool_index = 0;
				do {
					if(isFastQ){
						memcpy(info_list[pool_index].title, line, strlen(line) +1);
//						infile.getline(info_list[pool_index].title,dbfsequenceboundary);
						infile.getline(line,dbfsequenceboundary);
						lens_p[pool_index] = (int) strlen(line);
						lens_s[pool_index] = (int) strlen(line);
						memcpy( &lines_p[pool_index * dbfsequenceboundary], line, lens_p[pool_index] + 1);
						memcpy( &lines_s[pool_index * dbfsequenceboundary], line, lens_p[pool_index] + 1);
						infile.getline(info_list[pool_index].plus_head,dbfsequenceboundary);
						infile.getline(info_list[pool_index].plus_content,dbfsequenceboundary);
						if (!infile.eof())
						{
							infile.getline(line,dbfsequenceboundary);
						}

						pool_index ++;
						count++ ;

					}else {

						memcpy(temp_title, line, strlen(line) +1);
						while (!infile.eof() ){

							infile.getline(line,dbfsequenceboundary);
							count++;
							if (line[0] == '>') {
								break;
							}
							memcpy(info_list[pool_index].title, temp_title, strlen(temp_title) +1);
							lens_p[pool_index] = (int) strlen(line);
							lens_s[pool_index] = (int) strlen(line);
							memcpy( &lines_p[pool_index * dbfsequenceboundary], line, lens_p[pool_index] + 1);
							memcpy( &lines_s[pool_index * dbfsequenceboundary], line, lens_p[pool_index] + 1);
							pool_index ++;

							if (pool_index >= POOL_SIZE) {
								memcpy( line, temp_title, strlen(temp_title) +1);
								break;
							}
						}

					}
				 // count number of nucleotide sequences
				} while ( !infile.eof() && pool_index <POOL_SIZE);


				// Copy data from host to device;
				cudaMemcpy(d_lines, lines_p, dbfsequenceboundary * POOL_SIZE * sizeof(char), cudaMemcpyHostToDevice );
				cudaMemcpy(d_lens, lens_p, POOL_SIZE * sizeof(int), cudaMemcpyHostToDevice);

				// Launch kernel to process for POOL_SIZE of sequences lines_p
				kernel_findp<<< (pool_index -1)/BLOCK +1 , BLOCK>>>(d_lines, pool_index , dbfsequenceboundary , d_type, d_lens,
						 ptag[0] , mismatch, scangap , tail_length,  head_gap, d_true_identity, d_pos_list, d_pos_count);

				// Copy result back from device to host
				cudaMemcpy(lines_p, d_lines, dbfsequenceboundary * POOL_SIZE * sizeof(char), cudaMemcpyDeviceToHost );
				cudaMemcpy(true_identity_p, d_true_identity, POOL_SIZE * sizeof(float), cudaMemcpyDeviceToHost);
				cudaMemcpy(type_p, d_type, POOL_SIZE * sizeof(int), cudaMemcpyDeviceToHost);
				cudaMemcpy(pos_list_p, d_pos_list, 2 * POOL_SIZE * sizeof(position), cudaMemcpyDeviceToHost);
				cudaMemcpy(pos_count_p, d_pos_count, POOL_SIZE * sizeof(int), cudaMemcpyDeviceToHost);


				// Copy data from host to device;
				cudaMemcpy(d_lines, lines_s, dbfsequenceboundary * POOL_SIZE * sizeof(char), cudaMemcpyHostToDevice );
				cudaMemcpy(d_lens, lens_s, POOL_SIZE * sizeof(int), cudaMemcpyHostToDevice);

				// Launch kernel to process for POOL_SIZE of sequences lines_p
				kernel_finds<<< (pool_index -1)/BLOCK +1 , BLOCK>>>(d_lines, pool_index , dbfsequenceboundary , d_type, d_lens,
						 stag[0] , mismatch, scangap , tail_length,  head_gap, d_true_identity, d_pos_list, d_pos_count);

				// Copy result back from device to host
				cudaMemcpy(lines_s, d_lines, dbfsequenceboundary * POOL_SIZE * sizeof(char), cudaMemcpyDeviceToHost );
				cudaMemcpy(true_identity_s, d_true_identity, POOL_SIZE * sizeof(float), cudaMemcpyDeviceToHost);
				cudaMemcpy(type_s, d_type, POOL_SIZE * sizeof(int), cudaMemcpyDeviceToHost);
				cudaMemcpy(pos_list_s, d_pos_list, 2 * POOL_SIZE * sizeof(position), cudaMemcpyDeviceToHost);
				cudaMemcpy(pos_count_s, d_pos_count, POOL_SIZE * sizeof(int), cudaMemcpyDeviceToHost);

				// print result
				for (int lIdx = 0 ; lIdx < pool_index ; lIdx++) {
					char * lp = & lines_p[lIdx * dbfsequenceboundary];
					char * ls = & lines_s[lIdx * dbfsequenceboundary];

					if(output_s=="tl" && type_p[lIdx] == TYPE_0TR_1TL)
					{
						int temp_start_t = pos_list_p[lIdx].start;
						int temp_end_t = pos_list_p[lIdx].end;
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						if(detail==1)
						{
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<","<<"length="<<temp_end_t-temp_start_t+1<<","<<"identity="<<true_identity_p[lIdx]<<"%"<<")"<<endl;
						}
						else
						{
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<")"<<endl;
						}

						cout<< lp <<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}
					////////end tl outputs
					if(output_s=="tr" && type_p[lIdx] == TYPE_1TR_0TL)
					{
						int temp_start_t = pos_list_p[lIdx].start;
						int temp_end_t = pos_list_p[lIdx].end;
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						if(detail==1)
						{
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<","<<"length="<<temp_end_t-temp_start_t+1<<","<<"identity="<<true_identity_p[lIdx]<<"%"<<")"<<endl;
						}
						else
						{
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<")"<<endl;
						}

						cout<< lp <<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}

					//end tr
					if(output_s=="td" &&  type_p[lIdx] == TYPE_1TR_1TL)
					{
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						int temp_start_t1 = pos_list_p[lIdx].start;
						int temp_end_t1 = pos_list_p[lIdx].end;
						int temp_start_t2 = pos_list_p[lIdx + POOL_SIZE].start;
						int temp_end_t2 = pos_list_p[lIdx + POOL_SIZE].end;


						if(detail==1)
						{
							if (pos_count_p[lIdx] > 1)
								cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<","<<"length="<<temp_end_t1-temp_start_t1+1<<","<<"identity="<<true_identity_p[lIdx]<<"%"<<")"<<"("<<"type_p="<<ptag<<","<<"start="<<temp_start_t2+1<<","<<"end="<<temp_end_t2+1<<","<<"length="<<temp_end_t2-temp_start_t2+1<<","<<"identity="<<true_identity_p[lIdx+POOL_SIZE]<<"%"<<")"<<endl;

								 cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<","<<"length="<<temp_end_t1-temp_start_t1+1<<","<<"identity="<<true_identity_p[lIdx]<<"%"<<")"<<endl;

						}
						else
						{
							if (pos_count_p[lIdx] > 1)
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<")"<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t2+1<<","<<"end="<<temp_end_t2+1<<")"<<endl;
							else if (pos_count_p[lIdx] == 1)
								 cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<")"<<endl;
						}


						//change capital character
						cout<< lp <<endl;

						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}

					}
					//End td

					if(output_s=="nt" && type_p[lIdx]==TYPE_NONE)
					{

						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;
						}
						cout<<endl<<lp<<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}

					//////////////////////////////////////////////////////
					if(output_s=="al" && type_s[lIdx] == TYPE_0TR_1TL)
					{
						int temp_start_t = pos_list_s[lIdx].start;
						int temp_end_t = pos_list_s[lIdx].end;
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						if(detail==1)
						{
							cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<","<<"length="<<temp_end_t-temp_start_t+1<<","<<"identity="<<true_identity_s[lIdx]<<"%"<<")"<<endl;
						}
						else
						{
							cout<<" "<<"("<<"type="<<ptag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<")"<<endl;
						}

						cout<< ls <<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}
					////////end tl outputs
					if(output_s=="ar" && type_s[lIdx] == TYPE_1TR_0TL)
					{
						int temp_start_t = pos_list_s[lIdx].start;
						int temp_end_t = pos_list_s[lIdx].end;
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						if(detail==1)
						{
							cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<","<<"length="<<temp_end_t-temp_start_t+1<<","<<"identity="<<true_identity_s[lIdx]<<"%"<<")"<<endl;
						}
						else
						{
							cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t+1<<","<<"end="<<temp_end_t+1<<")"<<endl;
						}

						cout<< ls <<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}

					//end ar
					if(output_s=="ad" &&  type_s[lIdx] == TYPE_1TR_1TL)
					{
						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;

						}
						int temp_start_t1 = pos_list_s[lIdx].start;
						int temp_end_t1 = pos_list_s[lIdx].end;
						int temp_start_t2 = pos_list_s[lIdx + POOL_SIZE].start;
						int temp_end_t2 = pos_list_s[lIdx + POOL_SIZE].end;


						if(detail==1)
						{
							if (pos_count_s[lIdx] > 1)
								cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<","<<"length="<<temp_end_t1-temp_start_t1+1<<","<<"identity="<<true_identity_s[lIdx]<<"%"<<")"<<"("<<"type_s="<<stag<<","<<"start="<<temp_start_t2+1<<","<<"end="<<temp_end_t2+1<<","<<"length="<<temp_end_t2-temp_start_t2+1<<","<<"identity="<<true_identity_s[lIdx+POOL_SIZE]<<"%"<<")"<<endl;

							cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<","<<"length="<<temp_end_t1-temp_start_t1+1<<","<<"identity="<<true_identity_s[lIdx]<<"%"<<")"<<endl;

						}
						else
						{
							if (pos_count_s[lIdx] > 1)
								cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<")"<<"("<<"type_s="<<stag<<","<<"start="<<temp_start_t2+1<<","<<"end="<<temp_end_t2+1<<")"<<endl;
							else if (pos_count_s[lIdx] == 1)
								cout<<" "<<"("<<"type="<<stag<<","<<"start="<<temp_start_t1+1<<","<<"end="<<temp_end_t1+1<<")"<<endl;
						}


						//change capital character
						cout<< ls <<endl;

						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}

					}
					//End ad

					if(output_s=="na" && type_s[lIdx]==TYPE_NONE)
					{

						if(output_format == "fastq")
						{
							cout<< info_list[lIdx].title ;
						}
						if(output_format == "fasta")
						{
							info_list[lIdx].title[0] = '>' ;;
							cout<< info_list[lIdx].title ;
						}
						cout<<endl<<ls<<endl;
						//output quality line
						if(output_format == "fastq")
						{
							cout << info_list[lIdx].plus_head<<endl;
							cout << info_list[lIdx].plus_content<<endl;
						}
					}


					if(type_p[lIdx] == TYPE_0TR_1TL && type_s[lIdx]== TYPE_NONE)
					{
						t_left++;

					}
					if(type_p[lIdx] == TYPE_1TR_0TL  && type_s[lIdx]== TYPE_NONE)
					{
						t_right++;
					}
					if(type_p[lIdx]== TYPE_NONE && type_s[lIdx] == TYPE_0TR_1TL)
					{
						a_left++;
					}
					if(type_p[lIdx]== TYPE_NONE && type_s[lIdx] == TYPE_1TR_0TL)
					{
						a_right++;
					}
					if(type_p[lIdx] == TYPE_1TR_1TL && type_s[lIdx]== TYPE_NONE)
					{
						t_double++;
					}
					if(type_p[lIdx]== TYPE_NONE && type_s[lIdx] == TYPE_1TR_1TL)
					{
						a_double++;
					}
					if(type_p[lIdx] == TYPE_0TR_1TL && type_s[lIdx] == TYPE_1TR_0TL)
					{
						tl_ar++;
					}
					if(type_p[lIdx] == TYPE_1TR_0TL && type_s[lIdx] == TYPE_0TR_1TL)
					{
						al_tr++;
					}


				}


			}




			if(output_s=="n")
			{
				cout<<"[Total Reads]:"<<count/div<<endl;
				cout<<"[Total Reads with desired primary homopolymers](case1+case2-case5):"<<t_left + t_right - t_double<<endl;
				cout<<"[Percentage of primary homopolymers detection]:"<<((t_left + t_right - t_double)/(count/div))*100<<"%"<<endl;
				cout<<"[Total Reads with desired secondary homopolymers](case3+case4-case6):"<<a_left + a_right - a_double<<endl;
				cout<<"[Percentage of secondary homopolymers detection]:"<<((a_left + a_right - a_double)/(count/div))*100<<"%"<<endl;
				cout<<"case1:["<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<".................]:"<<t_left<<"\t"<<(t_left/(count/div))*100<<"%"<<endl;
				cout<<"case2:[................."<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<"]:"<<t_right<<"\t"<<(t_right/(count/div))*100<<"%"<<endl;
				cout<<"case3:["<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<".................]:"<<a_left<<"\t"<<(a_left/(count/div))*100<<"%"<<endl;
				cout<<"case4:[................."<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<"]:"<<a_right<<"\t"<<(a_right/(count/div))*100<<"%"<<endl;
				cout<<"case5:["<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<"........."<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<"]:"<<t_double<<"\t"<<(t_double/(count/div))*100<<"%"<<endl;
				cout<<"case6:["<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<"........."<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<"]:"<<a_double<<"\t"<<(a_double/(count/div))*100<<"%"<<endl;
				cout<<"case7:["<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<"........."<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<"]:"<<tl_ar<<"\t"<<(tl_ar/(count/div))*100<<"%"<<endl;
				cout<<"case8:["<<stag[0]<<stag[0]<<stag[0]<<stag[0]<<"........."<<ptag[0]<<ptag[0]<<ptag[0]<<ptag[0]<<"]:"<<al_tr<<"\t"<<(al_tr/(count/div))*100<<"%"<<endl;
			}

			//Free host memory
			delete [] info_list;
			delete [] lines_p;
			delete [] lines_s;
			delete [] temp_title;
			delete [] lens_p;
			delete [] lens_s;
			delete [] true_identity_p;
			delete [] true_identity_s;
			delete [] type_p;
			delete [] type_s;
			delete [] pos_list_p;
			delete [] pos_list_s;
			delete [] pos_count_p;
			delete [] pos_count_s;
			/// Free CUDA memory;
			cudaFree(d_pos_count);
			cudaFree(d_type);
			cudaFree(d_lens);
			cudaFree(d_lines);
			cudaFree(d_true_identity);
			cudaFree(d_pos_list);


		}
		infile.close();
	}

	return 0;
}

__global__ void kernel_findp(char * lines, int lines_count , int line_lenght, int * type,int * lens, char ptag ,int mismatch,
		int scangap ,int tail_length, int head_gap, float * true_identity, position * pos_list ,int * d_pos_count)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	if ( idx < lines_count) {
		char * line = &lines[idx * line_lenght];
		int len = lens[idx];
		int ap_num_T=0;
		int temp_start_t =0;
		int temp_end_t =0;
		double tl_num=0;
		double tr_num=0;
		int pos_count = 0;

		int temp_position_p=0;

		for(int i=0;i<len;i++)
		{
			if(line[i]==ptag && line[i+1]==ptag)
			{
				if(i-temp_position_p>mismatch+5)
				{
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;
				}
				int temp_mismatch=0;
				for(int k=0;k<=scangap+1;k++)
				{
					if(i+2+k>=len)
					{
						break;
					}
					if(line[i+2+k]!=ptag)
					{
						temp_mismatch++;
					}
				}
				if(temp_mismatch>mismatch  && (ap_num_T+2)>=(tail_length) )
				{
					ap_num_T=ap_num_T+2;
				}
				if(temp_start_t!=0 && temp_mismatch==0)
				{
					temp_end_t = i+1;
					ap_num_T++;
				}
				if(temp_start_t==0 && temp_mismatch==0)
				{
					temp_start_t=i;
					temp_end_t = i+1;
					ap_num_T++;
				}
				if(temp_start_t==0 && temp_mismatch<=mismatch)
				{

					temp_start_t=i;
					temp_end_t = i+1+temp_mismatch;
					ap_num_T=ap_num_T+temp_mismatch;
				}

				if(temp_start_t!=0 && temp_mismatch<=mismatch)
				{
					temp_end_t = i+1+temp_mismatch;
					ap_num_T=ap_num_T+temp_mismatch;
				}

				if(i+3>=len && temp_mismatch==0 && temp_start_t!=0 && ((ap_num_T+1)>=(tail_length+1)))
				{
					temp_end_t = i+2;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					ap_num_T++;
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;

						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;
						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;
				}
				if(i+3>=len && temp_mismatch<=mismatch && temp_mismatch>0 && temp_start_t!=0 && ((ap_num_T+1)>=(tail_length+1)))
				{
					temp_end_t = i+2;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					ap_num_T++;
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;
						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;

						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;

				}
				if((temp_mismatch>mismatch || (i-temp_position_p>mismatch+10)) && ((ap_num_T+1)>=(tail_length)))
				{
					if(temp_start_t==1 && line[0]==ptag)
					{
						temp_start_t=0;
					}
					temp_end_t = i+2;
					if(i-temp_position_p>mismatch+2)
					{
						temp_end_t=temp_position_p;
					}
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					if(line[temp_end_t]==ptag && line[temp_end_t+1]==ptag)
					{
						temp_end_t=i+3;
					}
					ap_num_T++;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i;
					}

					//fix the length
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;
						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;
						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}

					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;

				}
				if(temp_mismatch>mismatch && ((ap_num_T+1)<(tail_length)))
				{
					temp_start_t=0;type[idx] =
					temp_end_t = 0;
					ap_num_T=0;

				}
				temp_position_p=i;
			}
		}
		type[idx] = -1;
		if (tr_num == 0 && tl_num ==0)
			type[idx] = TYPE_NONE;
		if ( (tr_num == 0 && tl_num ==1) || (tr_num == 1 && tl_num ==0) ){
			d_pos_count[idx] = pos_count;
			if (tr_num == 0 && tl_num ==1)
				type[idx] = TYPE_0TR_1TL;
			else
				type[idx] = TYPE_1TR_0TL;

			int miss_match_count = 0 ;
			for (int i = pos_list[idx ].start ; i <= pos_list[idx].end; i++) {
				if (line[i] != ptag) {
					miss_match_count ++;
				}
				if (line[i] < 97)
					line[i]+= 32;
				else line[i]-=32;
			}
			true_identity[idx]=(1-miss_match_count/float(pos_list[idx].end-pos_list[idx].start+1))*100;
		}
		if ( tr_num == 1 && tl_num ==1 ) {
			d_pos_count[idx] = pos_count;
			type[idx] = TYPE_1TR_1TL;
			int miss_match_count = 0 ;
			for (int i = pos_list[idx ].start  ; i <= pos_list[idx].end; i++) {
				if (line[i] != ptag) {
					miss_match_count ++;
				}
				if (line[i] < 97)
					line[i]+= 32;
				else line[i]-=32;
			}
			true_identity[idx]=(1-miss_match_count/float(pos_list[idx].end-pos_list[idx].start+1))*100;
			if (pos_count > 1) {
				int miss_match_count2 = 0 ;
				for (int i = pos_list[idx + POOL_SIZE].start ; i <= pos_list[idx + POOL_SIZE].end; i++) {
					if (line[i] != ptag) {
						miss_match_count2 ++;
					}
					if (line[i] < 97)
						line[i]+= 32;
					else line[i]-=32;
				}
				true_identity[idx + POOL_SIZE]=(1-miss_match_count2/float(pos_list[idx + POOL_SIZE].end-pos_list[idx + POOL_SIZE].start+1))*100;
			}
		}
	}
}

__global__ void kernel_finds(char * lines, int lines_count , int line_lenght, int * type,int * lens, char ptag ,int mismatch,
		int scangap ,int tail_length, int head_gap, float * true_identity, position * pos_list ,int * d_pos_count)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	if ( idx < lines_count) {
		char * line = &lines[idx * line_lenght];
		int len = lens[idx];
		int ap_num_T=0;
		int temp_start_t =0;
		int temp_end_t =0;
		double tl_num=0;
		double tr_num=0;
		int pos_count = 0;

		int temp_position_p=0;

		for(int i=0;i<len;i++)
		{
			if(line[i]==ptag && line[i+1]==ptag)
			{
				if(i-temp_position_p>mismatch+5)
				{
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;
				}
				int temp_mismatch=0;
				for(int k=0;k<=scangap+2;k++)
				{
					if(i+2+k>=len)
					{
						break;
					}
					if(line[i+2+k]!=ptag)
					{
						temp_mismatch++;
					}
				}
				if(temp_mismatch>mismatch  && (ap_num_T+2)>=(tail_length) )
				{
					ap_num_T=ap_num_T+2;
				}
				if(temp_start_t!=0 && temp_mismatch==0)
				{
					temp_end_t = i+1;
					ap_num_T++;
				}
				if(temp_start_t==0 && temp_mismatch==0)
				{
					temp_start_t=i;
					temp_end_t = i+1;
					ap_num_T++;
				}
				if(temp_start_t==0 && temp_mismatch<=mismatch)
				{

					temp_start_t=i;
					temp_end_t = i+1+temp_mismatch;
					ap_num_T=ap_num_T+temp_mismatch;
				}

				if(temp_start_t!=0 && temp_mismatch<=mismatch)
				{
					temp_end_t = i+1+temp_mismatch;
					ap_num_T=ap_num_T+temp_mismatch;
				}

				if(i+3>=len && temp_mismatch==0 && temp_start_t!=0 && ((ap_num_T+1)>=(tail_length+1)))
				{
					temp_end_t = i+2;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					ap_num_T++;
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;

						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;
						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;
				}
				if(i+3>=len && temp_mismatch<=mismatch && temp_mismatch>0 && temp_start_t!=0 && ((ap_num_T+1)>=(tail_length+1)))
				{
					temp_end_t = i+2;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					ap_num_T++;
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;
						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;

						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}
					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;

				}
				if((temp_mismatch>mismatch || (i-temp_position_p>mismatch+12)) && ((ap_num_T+1)>=(tail_length)))
				{
					if(temp_start_t==1 && line[0]==ptag)
					{
						temp_start_t=0;
					}
					temp_end_t = i+2;
					if(i-temp_position_p>mismatch+1)
					{
						temp_end_t=temp_position_p;
					}
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i+1;
					}
					if(line[temp_end_t]==ptag && line[temp_end_t+1]==ptag)
					{
						temp_end_t=i+3;
					}
					ap_num_T++;
					if(line[temp_end_t]!=ptag)
					{
						temp_end_t=i;
					}

					//fix the length
					if(temp_end_t-temp_start_t+1>=tail_length)
					{
						if(temp_start_t<head_gap)
						{
							tl_num++;
						}
						if(temp_end_t>(len-head_gap))
						{
							tr_num++;
						}
						pos_list[idx +pos_count * POOL_SIZE].start = temp_start_t;
						pos_list[idx +pos_count * POOL_SIZE].end = temp_end_t;
						pos_count ++;
					}

					temp_start_t=0;
					temp_end_t = 0;
					ap_num_T=0;

				}
				if(temp_mismatch>mismatch && ((ap_num_T+1)<(tail_length)))
				{
					temp_start_t=0;type[idx] =
					temp_end_t = 0;
					ap_num_T=0;

				}
				temp_position_p=i + 1;
			}
		}
		type[idx] = -1;
		if (tr_num == 0 && tl_num ==0)
			type[idx] = TYPE_NONE;
		if ( (tr_num == 0 && tl_num ==1) || (tr_num == 1 && tl_num ==0) ){
			d_pos_count[idx] = pos_count;
			if (tr_num == 0 && tl_num ==1)
				type[idx] = TYPE_0TR_1TL;
			else
				type[idx] = TYPE_1TR_0TL;

			int miss_match_count = 0 ;
			for (int i = pos_list[idx ].start ; i <= pos_list[idx].end; i++) {
				if (line[i] != ptag) {
					miss_match_count ++;
				}
				if (line[i] < 97)
					line[i]+= 32;
				else line[i]-=32;
			}
			true_identity[idx]=(1-miss_match_count/float(pos_list[idx].end-pos_list[idx].start+1))*100;
		}
		if ( tr_num == 1 && tl_num ==1 ) {
			d_pos_count[idx] = pos_count;
			type[idx] = TYPE_1TR_1TL;
			int miss_match_count = 0 ;
			for (int i = pos_list[idx ].start  ; i <= pos_list[idx].end; i++) {
				if (line[i] != ptag) {
					miss_match_count ++;
				}
				if (line[i] < 97)
					line[i]+= 32;
				else line[i]-=32;
			}
			true_identity[idx]=(1-miss_match_count/float(pos_list[idx].end-pos_list[idx].start+1))*100;
			if (pos_count > 1) {
				int miss_match_count2 = 0 ;
				for (int i = pos_list[idx + POOL_SIZE].start ; i <= pos_list[idx + POOL_SIZE].end; i++) {
					if (line[i] != ptag) {
						miss_match_count2 ++;
					}
					if (line[i] < 97)
						line[i]+= 32;
					else line[i]-=32;
				}
				true_identity[idx + POOL_SIZE]=(1-miss_match_count2/float(pos_list[idx + POOL_SIZE].end-pos_list[idx + POOL_SIZE].start+1))*100;
			}
		}
	}
}

